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The 3dFIAT code simulates pyrolysis, ablation, and shape change of thermal protection 
materials and systems in three dimensions. The governing equations, which include energy 
conservation, a three-component decomposition model, and a surface energy balance, are 
solved with a moving grid system to simulate the shape change due to surface recession. This 
work is the first part of a code validation study for new capabilities that were added to 
3dFIAT. These expanded capabilities include a multi-block moving grid system and an 
orthotropic thermal conductivity model. This paper focuses on conditions with minimal 
shape change in which the fluid/solid coupling is not necessary. Two groups of test cases of 
3dFIAT analyses of Phenolic Impregnated Carbon Ablator in an arc -jet are presented. In 
the first group, axisymmetric iso-q shaped models are studied to check the accuracy of three- 
dimensional multi-block grid system. In the second group, similar models with various 
through-the-thickness conductivity directions are examined. In this group, the material 
thermal response is three-dimensional, because of the carbon fiber orientation. Predictions 
from 3dFIAT are presented and compared with arcjet test data. The 3dFIAT predictions 
agree very well with thermocouple data for both groups of test cases. 

Nomenclature 

A = area, m 2 

B' = ml p e u e C M , dimensionless mass blowing rate 

B a = pre-exponential constant, s" 1 

C H ,C M - Stanton numbers for heat and mass transfer 

c p = specific heat, J/kg-K 

E a = activation energy, J/kmol 

= outward pyrolysis mass flux, kg/m 2 -s 
= recovery enthalpy, J/kg 
= enthalpy, J/kg 

= partial heat of charring, defined in Eq.(4), J/kg 
= thermal conductivity, W/m-K 
= mass flux, kg/m 2 -s 
= radiative heat flux at surface, W/m 2 

= gas constant in Eq.(6), J/kmol-K; or residual in Eqs.(10)-(12) 

= rotation matrix 
= temperature, K 
t = time, s 

v = local grid velocity, m/s 

x,y,z = Cartesian coordinate system, m 

x\y\z' = Cartesian coordinate system in principal directions, m 
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= diffusion-coefficient-weighted species fraction 
= surface absorptance 
= rotation angles about x\y\z' directions 
= surface emissivity 
= volume fraction of resin 
= surface normal direction 
= blowing reduction parameter 
= total density, kg/m 3 

= original density of pyrolysis gas component, kg/m 3 
= residual density of pyrolysis gas component, kg/m 3 
= Stefan-Boltzmann constant, W/m 2 -K 4 
= mass fraction of virgin material, defined in Eq.(3) 

= decomposition reaction order 

= char 

= boundary-layer edge 
= pyrolysis gas 

= density component (A, B, and C) 

= surface species 
= virgin 
= wall 

= location of TPS inner surface 


I. Introduction 

The authors have developed a family of programs for analysis of ablative thermal protection system (TPS) 
materials. Fully Implicit Ablation and Thermal response code (FIAT), Two-dimensional Implicit Thermal response 
and Ablation code (TITAN), and 3 -dimensional Finite- volume alternatively directional Implicit Ablation and 
Thermal response code (3dFIAT) simulate the internal heat conduction, in-depth thermal decomposition, quasi- 
steady pyrolysis gas flow, and surface ablation of TPS materials in one, two, and three dimensions, respectively. 1-3 
FIAT is widely used by NASA and industry as the one-dimensional analysis and sizing tool for spacecraft TPS 
materials. TITAN can analyze problems with two-dimensional or axisymmetric geometry. In some cases, a two- 
dimensional analysis is inadequate, and a three-dimensional ablation code is required to perform a high fidelity 
simulation. The 3dFIAT program can analyze the thermal response of the entire heatshield of a space vehicle. The 
prediction of ablative heatshield response for a spacecraft entering the atmosphere with an angle of attack is such a 
case. The first version of 3dFIAT was a single -block finite- volume program, and it was not sufficiently flexible to 
model a complicated TPS/structure with multiple materials. However, by integrating 3dFIAT and MARC, 4 a 
simulation system was developed. This system can predict surface recession, shape change, in-depth pyrolysis, and 
internal thermal response for a three-dimensional TPS/structure system under general hypersonic entry conditions. 

The purposes of this paper are to demonstrate and validate the new capabilities that were added to the 3dFIAT 
program specifically for analysis of TPS materials. These expanded capabilities include a multiple -block moving 
grid system and a model for orthotropic thermal conductivity. A multiple -block system allows 3dFIAT to perform 
computations for some TPS/structure systems with complicated geometry and multiple materials without integrating 
with MARC. The properties of many thermal protection materials are transverse isotropic, which is a subset of 
orthotropic. For transverse isotropic materials, depending on fiber orientation, there is a primary “through-the- 
thickness” (TTT) direction with thermal conductivity that is different (and typically lower) from the value in the 
perpendicular plane (call the IP direction). The fibers are in the IP direction. Phenolic Impregnated Carbon Ablator 
(PICA) is one example of an orthotropic material. 5 

This is the first part of 3dFIAT validation effort. The results presented in this paper focus on conditions with 
surface recession but without significant shape change, such that the fluid/solid shape change coupling is not 
necessary. Conditions that require fluid/solid coupling will be discussed in a separate paper. Here, two groups of test 
cases of 3dFIAT analyses of PICA in an arc -jet are considered. In the first group, an axisymmetric iso-q shaped 
model is studied to check the accuracy of the three-dimensional multi-block moving grid system. In the second 
group, similar models with various TTT conductivity directions are examined. In this group, the material thermal 
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response is three-dimensional even if the surface geometry is symmetric, because of the carbon fiber orientation. For 
both groups, the predicted in-depth temperature histories are compared with available thermocouple data. 

II. Governing Equations 

The internal energy balance 3 * is a transient thermal conduction equation with additional pyrolysis terms: 


P c p 


dr 

dt 


V -(kVT)-(h g -h)V -m g - m g - Vh g + pc p v-VT 


( 1 ) 


The individual terms in Eq.(l) are interpreted as follows: rate of storage of sensible energy, net rate of thermal 
conductive heat flux, pyrolysis energy-consumption rate, net rate of energy convected by pyrolysis, and convection 
rate of sensible energy due to coordinate system movement. 

The specific heat is input as a function of temperature for both virgin and fully-charred material. In partially 
pyrolyzed zones (p c < p< p v ), the specific heat is obtained from the mixing rule 

C p =TC pv +(l-T) Cpc 


where 


Pv-Pc P 


( 3 ) 


The weighting variable r is the mass fraction of virgin material, in a hypothetical mixture of virgin material and 
char, which yields the correct local density. The thermal conductivity, k , is weighted in a similar manner as c p . 

The pyrolysis gas enthalpy, h g , is input as a function of temperature and pressure. The quantity h in Eq.(l), as 
defined below, is a function of temperature and is calculated from other input quantities. 

— _ p v h v ~p c h c 

P v ~Pc ( 4 ) 

Thermal and mechanical properties, including c p and k of virgin and char of many spacecraft heatshield materials, 
are available in the TPSX material properties database, accessible through the Internet . 6 

For pyrolyzing TPS materials, a standard three -component decomposition model is used . 7 The resin filler 
consists of two components A and B, and the reinforcing material C is the third component. The instantaneous local 
density of the composite is given by 


p= r 0^+p 6 )+(i-r)/? c (5) 

where the parameter T is the volume fraction of resin and is an input quantity. The three components decompose 
independently following the relation 


dp, 

dt 



(Pi-Prif 


K Pa J 


+ v-VPi 


( 6 ) 


where p ri is the residual or terminal density of component i, and p oi is the original density of component i . 

The motion of pyrolysis gas is assumed to be one -dimensional (in the surface-normal direction rj) and quasi- 
steady. Thus the mass flow rate of pyrolysis gas at the surface is calculated as 
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m 


( 7 ) 



where 7j 0 is located at the inner surface of the TPS. 

III. Boundary Conditions 

Conditions at the ablating surface are determined by convective and radiative heating and by surface 
thermochemical interactions with boundary -layer gases. The surface energy balance equation is written in general 
convective transfer-coefficient form as follows: 7 

+ m c h c + thgh g + a„q m - ae w T - q cw =0 

The first term in Eq. (8) represents the sensible convective heat flux, where h ew is the enthalpy of edge gas with 
frozen composition at the wall temperature. The sum of the second, third, and fourth terms in Eq. (8) is defined as 
the total chemical energy at the surface. The term containing Z* represents the transport of chemical energy 

associated with chemical reactions at the wall and in the boundary layer. The Z* driving forces for diffusive mass 
transfer include the effects of unequal diffusion coefficients. The fifth and sixth terms are the radiative heat fluxes 
absorbed and re-radiated by the wall, respectively, and the last term, q cw , represents the rate of heat conduction into 

the TPS. Here B f is the normalized mass blowing rate. The Aerotherm Chemical Equilibrium (ACE) 8 or Multi- 
component Ablation Thermochemistry (MAT) 9 codes can be used to generate tables of B' for charring materials. 

A blowing correction accounts for the reduction in transfer coefficients due to the transpiration of gases from 
pyrolysis and surface ablation into the boundary layer. The blowing rate correction equation for convective heat 
transfer is 


p e u e C H (H r -h ew ) + p e u e C N 


, * * \ T w 

2a z - z ) h - 

x je jw' j 


-B'h 


C H _ ln(l + 2Z B') 

~ 2 Jb' (9) 

where X is the blowing reduction parameter, C H is the heat transfer coefficient for the ablating surface, and C m is 

the heat transfer coefficient for the non-ablating surface. With A = 0.50, Eq. (9) reduces to the classical blowing 
correction for laminar flow. 10 

The time dependent temperature and density boundary conditions at the interface of two blocks are estimated 
using a linear interpolation routine based on the value at the center of boundary cells. 

IV. Orthotropic Thermal Conductivity 

For materials with orthotropic thermal conductivity, the thermal conductivities in three principle directions must 
be defined. The thermal conductivity of many thermal protection materials is transverse isotropic, which is a subset 
of orthotropic thermal conductivity. For transverse isotropic materials, there is a primary through-the-thickness 
(TTT) direction with thermal conductivity that is different (and often lower) than the value in the perpendicular 
plane (called the in-plane or IP direction). PICA is one example of a transverse isotropic material. To define a 
transverse isotropic material in 3dFIAT, the thermal conductivities in the principle orthogonal directions (TTT and 
IP) are input. For a general three-dimensional body rotation from the principle directions (x’,y’,z’) to the local 
Cartesian coordinates (x,y,z), the rotation matrix R may be expressed in terms of three angles (y, (3, and a), 
y, (3, and a are the counterclockwise rotations about x’, y’, and z’ axis, respectively. 
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R = R z .(a)R y .(/3)R x .(y) 


R z -(«) = 
R y(P) = 


cosa - sin a 0 
sin a cosa 0 

0 0 1 

cos /? 0 sin/? 

0 1 0 
-sin/? 0 cos/? 


R Ay) 


i o o 

0 cos y -sin y 
0 sin^ cos y 


( 10 ) 


and 


R = 


cosacos/? 
sinacos/? 
-sin /? 


cosasin /? sin y -sin a cos / 
sinasin /?sin 7 + cosacosf 
cos/?sinf 


cosasin /?cosf + sinasinf 
sinasin /?cosf- cosasin y 

COS/?COSf 


( 11 ) 
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K X y 
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0 

K IP 
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Kzy 

K Z z 


0 

0 

K ip _ 


( 12 ) 


For typical arcjet models, we may set both y and (3 equal to 0, and a is the angle in the x’-y’ plane between the TTT 
direction (x’) and the axis of symmetry (x). The rotation matrix, R, becomes a single rotation about the z-axis, R z >. 
The local Cartesian coordinates and the principal directions of an isotropic transverse material used for the 
simulation of isq-q shaped arc -jet models are shown in Fig. 1. For baseline PICA material, the TTT direction is 
parallel to the axis of symmetry, and a is equal to 0°. Thus K ^ = K TT , K yy = K zz = K IP , and the rest of elements 

in conductivity tensor are zero. 


V. Test Cases 

The computations presented in this paper for 3dFIAT validation focus on analysis of iso-q shaped models 
typically used in the arc -jet testing conducted at NASA Ames Research Center. It has been shown in the previous 
work that shape change due to surface recession on an iso-q model is insignificant after its exposure to arc -jet 
stream. 11 This is because the surface recession over the front surface of model is fairly uniform. Consequently, the 
heating and pressure profiles over the model surface remain unchanged during its exposure in the arc -jet stream. The 
material thermal response simulation of an iso-q shaped model can thus be performed without flow/solid shape 
change coupling. 

There are two groups of cases studied in this work. The first group has arc -jet models with baseline PICA, 
which has the TTT direction aligned with the flow direction (a = 0°). In this group, the material thermal response is 
axisymmetric. The predictions using a three-dimensional multi -block moving grid system are compared with the 
available thermocouple data to check the accuracy of the grid system. The second group has models with somewhat 
different thickness profiles and a equal to 0°, 45°, 70° and 90°. These cases are used to examine the orthotropic 
thermal conductivity model implemented in the code for three-dimensional thermal diffusion. The surface geometry 
of the test models is axisymmetric. However, the material thermal response is three-dimensional if a is non-zero. 
The predictions made in the second group are also compared with thermocouple data to study if the code is self- 
consistent and accurate. 

A. Group I 

The material map for the arc -jet model studied in the first group is shown in Fig. 2. The geometry is an iso-q- 
shaped PICA model with 10.16-cm diameter, and the model holder is made of LI-2200. The initial surface geometry 
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of test model is axisymmetric. Figure 3 is the computational grid system generated for Group I simulation. We use a 
three -block no-singular-line grid system. The grid size is 70x87x19 (block 1), 47x155x19 (block 2), and 6x132x11 
(block 3). Based on our experience, computations using an axisymmetric three-dimensional grid system with a 
singular line may result in anomalous predictions around the singular line. The purpose of the third block of the grid 
is to avoid the presence of a singular line along the axis of symmetry. 

Two test cases, Gla and Gib, are studied in the first group. The arc -jet exposure time, stagnation heat flux, and 
stagnation pressure for both cases are listed in Table I. The TTT direction (x’) of PICA material is parallel to the arc 
flow direction (x). This is called the baseline PICA configuration. The coordinate systems used in the baseline 
configuration are depicted in Fig. 4. The principle directions (x’, y’, z’) of thermal conductivity tensor coincide with 
the directions of local coordinates (x, y, z). Thus a is equal to 0°. The surface heating and pressure distributions, 
applied as the boundary conditions for 3dFIAT, are predicted by DPLR. 12 As discussed earlier, these conditions are 
assumed to be independent of time during the exposure to arc stream, because shape change is negligibly small. The 
predicted heating and pressure profiles over the model surface for case Gla are shown in Fig. 5. 


Table I: Arc -jet conditions 


Case 

Exposure Time (sec) 

Heat Flux (W/cm 2 ) 

Pressure (kPa) 

Gla 

42 

246 

8.5 

Gib 

60 

169 

5.0 

G2 

45 

768 

46.7 


The comparison of in-depth temperature history between thermocouple data and prediction by 3dFIAT for Case 
Gla is presented in Figs. 6a to 6c. The thermocouples, TCI to TC6, are located at the centerline of iso-q model 
(radius = 0) and at various depths. TC7, TC8, TC9, and TC10 are off the centerline (radius > 0). The location of 
each thermocouple is listed in Table II. “Depth” is the distance along the centerline from the initial stagnation point, 
and “Radius” is the distance to the centerline of model. The predictions (black lines) and the thermocouple data (red 
lines) are in excellent agreement for all the thermocouples. Similar comparisons for Case Gib are presented in Figs 
7a to 7c. Again, the agreement between computation and data is good. This agreement indicates that the three- 
dimensional multi-block no-singular-line grid system is correctly implemented in the 3dFIAT code. The stagnation 
point recession predicted by 3dFIAT is exactly the same as that by TITAN and FIAT, because the ablation models 
adopted in these three codes are identical. Thus, the comparison of predicted recession with data is not discussed 
here. The same computations performed by TITAN and FIAT can be found in Reference 5. Figure 8 shows the 
predicted temperature contours on the plane of symmetry (x-y plane) for Case Gla at time equal to 600 sec. As 
expected, the thermal response of the PICA model is axisymmetric. 


Table II: The locations of thermocouples for Group I 


Group I 

Depth (cm) 

Radius (cm) 

TCI 

0.38 

0 

TC2 

0.76 

0 

TC3 

1.14 

0 

TC4 

1.52 

0 

TC5 

2.29 

0 

TC6 

3.05 

0 

TC7 

2.29 

4.44 

TC8 

2.29 

3.81 

TC9 

2.29 

2.54 

TC10 

3.05 

4.44 


B. Group II 

The surface geometry of the arc -jet model for the second group is the same as that for the first group, but the 
thickness of the PICA material is different. The thickness of PICA along the centerline is 3.49 cm for this group, and 
it is 4.13 cm for the first group. The material map of Group II is shown in Fig. 9. The no -singular-line, three-block 
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computational grid system is presented in Fig. 10. The grid block sizes are the same as those for Group I. In this 
group, the iso-q shaped PICA models with various TTT angles (0°, 45°, 70°, and 90°) are studied. The arc -jet stream 
conditions are listed at row G2 in Table I. 

Figure 1 1 is the comparison of stagnation point total recession between 3dFIAT prediction and measured data. 
The predicted total recession (1.2 cm) is slightly lower than the measured value. The difference between prediction 
and data is less than 6%. The comparison between thermocouple data (TC11 to TCI 5) and 3dFIAT predictions at 
TTT directions of 0°, 45°, 70°, and 90° are presented in Figs. 12a to 12d, respectively. The locations of 
thermocouples are listed in Table III. Generally speaking, data and predictions are in good agreement for all the TTT 
directions. The only noticeable discrepancy between data and computation is at the location of TCI 5 for a TTT 
angle equal to 45°, in which the 3dFIAT prediction is higher than the data after 60 sec. The maximum difference 
reaches around 60° K at around 140 sec. Figure 13 shows the temperature history of TCI 5 for various TTT 
directions. TCI 5 is located at about 3.0 cm from the initial stagnation point. The computation did not catch the 
temperature rise that appears in the first 40 sec of the test. This early temperature rise has been seen in other carbon- 
phenolic materials. However, the cause of this phenomenon is still not clear. The predictions indicate that increasing 
the TTT angle should increase the temperature at the location of TCI 5. Predictions and data agree well except that 
for the TTT angle equal to 45°. It appears that the TCI 5 temperature reading for angle of 45° is in error, because its 
temperature reading ran below that of a = 0° after 140 sec. 

The temperature contours of the front surface and the plane of symmetry at time equal to 400 sec are shown in 
Figs. 14a to 14d. For TTT angles equal to 45° and 70°, the temperature at the lower half of the front surface clearly is 
higher than that at the upper half. This result is obtained because the carbon fiber orientation makes the upper corner 
cool down more effectively than the lower corner. At angle equal to 45°, the upper half front surface and the lower 
half front surface have the greatest temperature difference. For TTT angle of 90°, the TTT direction is parallel to the 
y axis, and thus the temperature contours of the lower half and the upper half are symmetric to the x-z plane 
(horizontal surface). The temperature distributions of the right half and the left half are also symmetric because of 
the geometry symmetry. However, the temperature contours are not symmetric to the axis of symmetry (x axis). 


Table III: The locations of thermocouples for Group II 


Group II, TTT angle 

0° 

45° 

o 

O 

90° 


Depth (cm) 




TC11 

0.59 

0.61 

0.63 

0.62 

TC12 

1.04 

0.97 

1.00 

1.00 

TC13 

1.41 

1.38 

1.37 

1.38 

TC14 

1.80 

1.74 

1.76 

1.75 

TC15 

3.04 

3.02 

3.03 

3.04 


VI. Conclusions 

A code validation study for new capabilities added to 3dFIAT was performed. These expanded capabilities 
include a multi-block moving grid system and an orthotropic thermal conductivity model. This paper focuses on 
conditions with surface recession but without significant shape change such that the fluid/solid coupling is not 
required. Two groups of test cases of 3dFIAT analyses for PICA model in an arc -jet stream were presented. In the 
first group, an axisymmetric iso-q shaped model was studied to check if three-dimensional multi-block moving grid 
system was properly implemented. A three-block no -singular-line grid system was used in this simulation. The 
agreement between 3dFIAT prediction and data was excellent for both the centerline and off-the-centerline 
thermocouples. In the second group, similar PICA models with various Through -The-Thickness (TTT) conductivity 
directions (0°, 45°, 70°, and 90°) were examined to check the implementation of orthotropic thermal conductivity 
model. Predictions from 3dFIAT were presented and also compared with arc -jet data. The difference on total surface 
recession between prediction and measurement is within 6%. The 3dFIAT predictions agreed well with the 
thermocouple data for all the TTT angles. It was found that the thermocouple temperature reading from TCI 5 for 
TTT angle equal to 45° was inconsistent with the readings from the rest of TTT angles, and was lower than the 
prediction by 3dFIAT. 
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Figure 1. Local coordinate system (x, y, z) and 
principle directions of conductivity (x% y’, z’). 


PICA LI-2200 



Figure 2. Material map for group I. 



Figure 3. Three-block grid system for group I. 



Figure 4. Coordinate systems for group I. 
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Figure 5. Surface heating and pressure profiles 
for Case IGa. 


Figure 6b. Predictions vs. data moving outward 
from model centerline at 2.29 cm depth, Case Gla. 
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Figure 6a. Predictions vs. data along model 
centerline, Case Gla. 



Figure 6c. Predictions vs. data moving outward 
from model centerline at 3.05 cm depth, Case Gla. 
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Figure 7a. Predictions vs. data along model 
centerline, Case Gib. 



Time (sec) 


Figure 7c. Predictions vs. data moving outward 
from model centerline at 3.04 cm depth, Case Gib. 




Figure 7b. Predictions vs. data moving outward Figure 8. Temperature contours at the plane of 

from model centerline at 2.29 cm depth, Case Gib. symmetry at time = 600 s. 
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Air Gap 


Figure 9, Material map for Group II. 



Figure 11. Stagnation point recession. 



Block 3 



Figure 10. Three-block system for Group II. Figure 12a. Prediction vs. data along model 

centerline; a = 0, Case G2. 
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Figure 12b. Prediction vs. data along model 
centerline; a = 45, Case G2. 



Figure 12c. Prediction vs. data along model 
centerline; a = 70, Case G2 



Fige 12d. Prediction vs. data along model 
centerline; a = 90, Case G2. 



Figure 13. Comparison of temperature history at 
TC15. 
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Figure 14a. Temperature contours at time = 400 s. 


Figure 14c. Temperature contours at time = 400 s. 



Figure 14b. Temperature contours at time = 400 s. Figure 14d. Temperature contours at time = 400 s. 
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